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METHOD AND APPARATUS FOR ADAPTIVE REAL-TIME 
SIGNAL CONDITIONING, PROCESSING, ANALYSIS, 
QUANTIFICATION, COMPARISON, AND CONTROL 

This application claims the benefit of United States Provisional Patent Application 
No.60/416,562 filed on October 7, 2002, which is incorporated herein by reference in its 
5 entirety. 

COPYRIGHT NOTIFICATION 

Portions of this patent application contain materials that are subject to copyright protection. 
The copyright owner has no objection to the facsimile reproduction by anyone of the patent 
10 document or the patent disclosure, as it appears in the Patent and Trademark Office patent 
file or records, but otherwise reserves all copyright rights whatsoever. 

TECHNICAL FIELD 

The present invention relates to methods, processes and apparatus for real-time measuring 
15 and analysis of variables. In particular, it relates to adaptive real-time signal conditioning, 
processing, analysis, quantification, comparison, and control. This invention also relates to 
generic measurement systems and processes, that is, the proposed measuring arrangements 
are not specially adapted for any specific variables, or to one particular environment. This 
invention also relates to methods and corresponding apparatus for measuring which extend 
20 to different applications and provide results other than instantaneous values of variables. The 
invention further relates to post-processing analysis of measured variables and to statistical 
analysis. 
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BACKGROUND ART 

Due to the rapid development of digital technology since the 1950's, the development of 
analog devices has been essentially squeezed out to the periphery of data acquisition equipment 
only. It could be argued that the conversion to digital technology is justified by the flexibility, 

5 universality, and low cost of modern integrated circuits. However, it usually comes at the price 
of high complexity of both hardware and software implementations. The added complexity 
of digital devices stems from the fact that all operations must be reduced to the elemental 
manipulation of binary quantities using primitive logic gates. Therefore, even such basic 
operations as integration and differentiation of functions require a very large number of such 

10 gates and/or sequential processing of discrete numbers representing the function sampled at 
many points. The necessity to perform a very large number of elemental operations limits 
the ability of digital systems to operate in real time and often leads to substantial dead time 
in the instruments. On the other hand, the same operations can be performed instantly in 
an analog device by passing the signal representing the function through a simple RC circuit. 

15 Further, all digital operations require external power input, while many operations in analog 
devices can be performed by passive components. Thus analog devices usually consume much 
less energy, and are more suitable for operation in autonomous conditions, such as mobile 
communication, space missions, prosthetic devices, etc. 

It is widely recognized (see, for example, Paul and Hiiper (1993)) that the main obstacle to 

20 robust and efficient analog systems often lies in the lack of appropriate analog definitions and 
the absence of differential equations corresponding to known digital operations. When proper 
definitions and differential equations are available, analog devices routinely outperform cor- 
responding digital systems, especially in nonlinear signal processing (Paul and Hiiper, 1993). 
However, there are many signal processing tasks for which digital algorithms are well known, 

25 but corresponding analog operations are hard to reproduce. One example, which is widely 
recognized to fall within this category, is related to the use of signal processing techniques 
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based on order statistics 1 . Order statistic (or rank) filters are gaining wide recognition for 
their ability to provide robust estimates of signal properties and are becoming the filters of 
choice for applications ranging from epileptic seizure detection (Osorio et aL, 1998) to image 
processing (Kim and Yaroslavsky, 1986). However, since such filters work by sorting, or order- 

5 ing, a set of measurements their implementation has been constrained to the digital domain. 
As pointed out by some authors (Paul and Huper, 1993, for example), the major problem in 
analog rank processing is the lack of an appropriate differential equation for 'analog sorting'. 
There have been several attempts to implement such sorting and to build continuous-time rank 
filters without using delay lines and/or clock circuits. Examples of these efforts include optical 

10 rank filters (Ochoa et aL, 1987), analog sorting networks (Paul and Huper, 1993; Opris, 1996), 
and analog rank selectors based on minimization of a non-linear objective function (Urahama 
and Nagao, 1995). However, the term 'analog' is often perceived as only 'continuous- time', 
and thus these efforts fall short of considering the threshold continuity, which is necessary for a 
truly analog representation of differential sorting operators. Even though Ferreira (2000, 2001) 

15 extensively discusses threshold distributions, these distributions are only piece wise-continuous 
and thus do not allow straightforward introduction of differential operations with respect to 
threshold. 

Nevertheless, fuelled by the need for robust filters that can operate in real time and on a low 
energy budget, analog implementation of traditionally digital operations has recently gained 

20 in popularity aided by the rapid progress in analog Very Large Scale Integration (VLSI) tech- 
nology (Mead, 1989; Murthy and Swamy, 1992; Kinget and Steyaert, 1997; Lee and Jen, 1993). 
However, current efforts to implement digital signal processing methods in analog devices still 
employ an essentially digital philosophy. That is, a continuous signal is passed through a 
delay line which samples the signal at discrete time intervals. Then the individual samples 

25 are processed by a cascade of analog devices that mimic elemental digital operations (Vlassis 
et aL, 2000). Such an approach fails to exploit the main strength of analog processing, which 
is the ability to perform complex operations in a single step without employing the 'divide 
and conquer' paradigm of the digital approach. 

1 See, for example, Arnold et al. (1992) for the definitions and theory of order statistics. 
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Perhaps the most common digital waveform device is the analog-to-digitaJ converter 
(ADC). Among the salient characteristics of ADCs are their sampling frequency, measure- 
ment resolution, power dissipation, and system complexity. Sampling frequency is typically 
dictated by the signal of interest and/or the requirements of the application. As the frequency 

5 content of the signal of interest increases and the sampling frequency increases, resolution 
decreases both in terms of the absolute number of bits available in an ADC and in terms of 
the effective number of bits (ENOB), or accuracy, of the measurement. Power needs typically 
increase with increasing sampling frequency. The system complexity is increased if continu- 
ous monitoring of an input signal is required (real-time operation). As an example, high-end 

10 oscilloscopes can capture fast transient events, but are limited by record length (the number 
of samples that can be acquired) and dead time (the time required to process, store, or dis- 
play the samples and then reset for more data acquisition). These limitations affect any data 
acquisition system in that, as the sampling frequency increases, resources will ultimately be 
limited at some point in the processing chain. In addition, the higher the acquisition speed, 

15 the more negative effects such as clock crosstalk, jitter, and synchronization issues combine to 
reduce system performance. 

It is highly desirable to extract signal characteristics or preprocess data prior to digitization 
so that the requirements on the ADCs are reduced and higher quality data can be obtained. 
In the past, one common technique was to use an analog memory to sample a fast signal and 

20 then the analog memory would be clocked out at a low speed and digitized with a moderately 
high resolution ADC. While this technique works, it suffers from significant degradation due 
to clock feedthrough, non-linear effects of the analog memories chosen, and limited record 
length. Another technique used is to multiplex a high-speed signal to a number of lower speed 
but higher resolution ADCs using an interleaved clock. Again, the technique works but never 

25 realizes the best performance of a single channel due to the high clock noise and inevitable 
differences in processing channels. 

The introduction of the Analysis of VAriabies Through Analog Representation (AVATAR) 
methodology (see Nikitin and Davidchack (2003a) and US Patent Application 09/921,524, 
which are incorporated herein by reference in their entirety) is aimed to address many aspects 
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of modern data acquisition and signal processing tasks by offering solutions that combine 
the benefits of both digital and analog technology, without the drawbacks of high cost, high 
complexity, high power consumption, and low reliability. The AVATAR methodology is based 
on the development of a new mathematical formalism, which takes into consideration the 

5 limited precision and inertial properties of real physical measurement systems. Using this 
formalism, many problems of signal analysis can be expressed in a content-sentient form 
suitable for analog implementation. Specific devices for a wide variety of signal processing tasks 
can be built from a few universal processing units. Thus, unlike traditional analog solutions, 
AVATAR offers a highly modular approach to system design. Most practical applications of 

10 AVATAR, however, are far from obvious, and their development requires technical solutions 
unavailable in the prior art. For example, AVATAR introduces the definitions of analog 
filters and selectors. Nonetheless, the practical implementations of these filters offered by 
AVATAR are often unstable and suffer from either lack of accuracy or lack of convergence 
speed, and thus are unsuitable for real-time processing of nonstationary signals. Another 

15 limitation of AVATAR lies in the definition of the threshold filter. Namely, a threshold filter 
in AVATAR depends only on the difference between the displacement and the input variables, 
and is expressed as a scalar function of only the displacement variable, which limits the 
scope of applicability of AVATAR. As another example, the analog counting in AVATAR is 
introduced through modulated density, and thus the instantaneous counting rate is expressed 

20 as a product of a rectified time derivative of the signal and the output of a probe. Even 
though this definition theoretically allows counting without dead time effect, its practical 
implementations are cumbersome and inefficient. 



5 



Atty. Ref. No. : Nikitin 03 01 
Inventor: Nikitin 
US Express: EK948840967US 



DISCLOSURE OF INVENTION 
BRIEF SUMMARY OF THE INVENTION 



The present invention, collectively designated as Adaptive Real-Time Embodiments for Mul- 
tivariate Investigation of Signals (ARTEMIS), overcomes the shortcomings of the prior art 

5 by directly processing the data in real-time in the analog domain, identifying the events of 
interest so that continuous digitization and digital processing is not required, performing 
direct, noise- resistant measurements of the salient signal characteristics, and outputting 
a signal proportional to these characteristics that can be digitized without the need for 
highspeed front-end sampling. 

10 In the face of the overwhelming popularity of digital technology, simple analog designs 
are often overlooked. Yet they often provide much cheaper, faster, and more efficient so- 
lutions in applications ranging from mobile communication and medical instrumentation to 
counting detectors for high-energy physics and space missions. The current invention, col- 
lectively designated as Adaptive Real-Time Embodiments for Multivariate Investigation of 

15 Signals (ARTEMIS), explores a new mathematical formalism for conducting adaptive content- 
sentient real-time signal processing, analysis, quantification, comparison, and control, and for 
detection, quantification, and prediction of changes in signals. The method proposed herein 
overcomes the limitations of the prior art by directly processing the data in real-time in the 
analog domain, identifying the events of interest so that continuous digitization and digital 

20 processing is not required, performing a direct measurement of the salient signal characteristics 
such as energy and arrival rate, and outputting a signal proportional to these characteristics 
that can be digitized without the need for highspeed front-end sampling. In addition, the 
analog systems can operate without clocks, which reduces the noise introduced into the data. 
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A simplified diagram illustrating multimodal analog real-time signal processing is shown in 
figure 1. The process comprises the step of Threshold Domain Filtering in combination with 
at least one of the following steps: (a) Multimodal Pulse Shaping, (b) Analog Rank Filtering, 
and (c) Analog Counting. 

5 Threshold Domain Filtering is used for separation of the features of interest in a signal from 
the rest of the signal. In terms of a threshold domain, a 'feature of interest' is either a point 
inside of the domain, or a point on the boundary of the domain. A typical Threshold Domain 
Filter can be composed of (asynchronous) comparators and switches, where the compara- 
tors operate on the differences between the components of the incoming variable(s) and the 

10 corresponding components of the control variable(s). For example, for the domain defined 
as a product of two ideal comparators represented by the Heaviside unit step function 6{x), 
V = 0[x(t) — D] 9[x(t)] (with two control levels, D and zero), a point inside (that is, V — 1) 
corresponds to a positive-slope signal of the magnitude greater than D, and the stationary 
points of x(t) above the threshold D can be associated with the points on the boundary of 

15 this domain. 

Multimodal Pulse Shaping can be used for embedding the incoming signal into a threshold 
space and thus enabling extraction of the features of interest by the Threshold Domain Filter- 
ing. A typical Multimodal Pulse Shaper transforms at least one component of the incoming 
signal into at least two components such that one of these two components is a (partial) 
20 derivative of the other. For example, for identification of the signal features associated with 
the stationary points of a signal x(t), the Multimodal Pulse Shaping is used to output both 
the signal x(t) and its time derivative x(t). 

Analog Rank Filtering can be used for establishing and maintaining the analog control lev- 
els of the Threshold Domain Filtering. It ensures the adaptivity of the Threshold Domain 
25 Filtering to changes in the measurement conditions, and thus the optimal separation of the 
features of interest from the rest of the signal. For example, the threshold level D in the 
domain V = 9{x(i) — D) 9[x(t)] can be established by means of Analog Rank Filtering to sep- 
arate the stationary points of interest from those caused by noise. Note that the Analog Rank 
Filtering outputs the control levels indicative of the salient properties of the input signal(s), 
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and thus can be used as a stand-alone embodiment of ARTEMIS for adaptive real-time signal 
conditioning, processing, analysis, quantification, comparison, and control, and for detection, 
quantification, and prediction of changes in signals. 

Analog Counting can be used for identification and quantification of the crossings of the 
5 threshold domain boundaries, and its output(s) can be either the instantaneous rate(s) of 
these crossings, or the rate(s) in moving window of time. A typical Analog Counter consists 
of a time-differentiator followed by a rectifier, and an optional time-integrator. 

Further scope of the applicability of the invention will be clarified through the detailed 
description given hereinafter. It should be understood, however, that the specific examples, 
10 while indicating preferred embodiments of the invention, axe presented for illustration only. 
Various changes and modifications within the spirit and scope of the invention should be- 
come apparent to those skilled in the art from this detailed description. Furthermore, all the 
mathematical expressions and the examples of hardware implementations are used only as a 
descriptive language to convey the inventive ideas clearly, and are not limitative of the claimed 
15 invention. 
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BRIEF DESCRIPTION OF FIGURES 

Figure 1. Simplified diagram of a multimodal analog system for real-time signal processing. 
Figure 2. Representative step (a) and impulse (b) responses of a continuous comparator. 
5 Figure 3. Defining output D q (t) of a rank filter as a level curve of the distribution function 
*(Z> f t). 

Figure 4. Amplitude and counting densities computed for the fragment of a signal shown at 
the top of the figure. 

Figure 5. Comparison of the rates measured with the boxcar (thin solid line) and the 
10 'triple-integrator' test function (n = 3 in equation (14), thick line) with r — T/6. The 
respective test functions are shown in the upper left corner of the figure. 
Figure 6. Illustration of bimodal pulse shaping used for the amplitude and timing measure- 
ments of a short-duration event. 

Figure 7. A simplified principal schematic of a signal processing system for a two-detector 
15 particle telescope. 

Figure 8. Example of analog measurement of instantaneous rate of signal's extrema. 
Figure 9. Example of analog counting of coincident maxima. 

Figure 10. Coincident counting in a time-of-flight window: (a). 'Ideal' instrument; (b) 'Re- 
alistic' instrument. 

20 Figure 11. Principle schematic of implementation of an adaptive real-time rank filter given 
by equation (24). 

Figure 12. Simplified diagram of an Adaptive Analog Rank Filter (AARF). 
Figure 13. Principle schematic of a 3-comparator implementation of AARF. 
Figure 14. Principle schematic of a delayed comparator. 
25 Figure 15. Principle schematic of an averaging comparator. 

Figure 16. Comparison of the performance of the analog rank filter given by equation (24) 
to that of the 'exact' quantile filter in a boxcar moving window of width T. 
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Figure 17. Example of using an internal reference (baseline) for separating signal from noise. 
Figure 18. Principle schematic of establishing baseline by means of quartile (trimean) filter. 
Figure 19. Principle schematic of implementation of an adaptive real-time rank selector 
given by equation (28). 

Figure 20. Simplified diagram of an Adaptive Analog Rank Selector (AARS). 

Figure 21. Removing static and dynamic impulse noise from a monochrome image by 

AARSs. 

Figure 22. Comparison of quartile outputs (q = 1/4, 1/2, and 3/4) of ARF and EARL for 
amplitudes and counts. 

Figure 23. Simplified principle schematic of a Bimodal Analog Sensor Interface System 
(BASIS). 

Figure 24. Simplified principle schematic of an Integrated Output Module of BASIS. 
Figure 25. Simulated example of the performance of BASIS used with a PMT. 
Figure 26. Modification of BASIS used for detection of onsets/offsets of light pulses. 
Figure 27. Principle schematic of a monoenergetic Poisson pulse generator. 
Figure 28 Simulated performance of a monoenergetic Poisson pulse generator. 
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TERMS AND DEFINITIONS 
WITH ILLUSTRATIVE EXAMPLES 

For convenience, the essential terms used in the subsequent detailed description of the 
invention are provided below, along with their definitions adopted for the purpose of this 
disclosure. Some examples clarifying and illustrating the meaning of the definitions are also 

5 provided. Note that the sections and equations in this part are numbered separately from the 
rest of the disclosure. Additional explanatory information on relevant terms and definitions 
can be found in US Patent Application 09/921,524 and US Provisional Patent Application 
No.60/416,562, which are incorporated herein by reference in their entirety. Some other terms 
and their definitions which might appear in this disclosure will be provided in the detailed 

10 description of the invention. 

D-l Continuous comparators and probes 

Consider a simple measurement process whereby a signal x(t) is compared to a threshold 
value D. The ideal measuring device would return '0' or T depending on whether x(t) is 

15 larger or smaller than D. The output of such a device is represented by the Heaviside unit 
step function 9 [D — x(t)], which is discontinuous at zero. However, the finite precision of real 
measurements inevitably introduces uncertainty in the output whenever x{t) « D. To describe 
this property of a real measuring device, we represent its output by a continuous function 
f AD [D — x(t)], where the width parameter AD characterizes the threshold interval over which 

20 the function changes from '0' to T and, therefore, reflects the measurement precision level. 
We call F& D (D) the threshold step response of a continuous comparator. Because of the 
continuity of this function, its derivative f± D {D) = dF AD /dD exists everywhere, and we 
call it the comparator's threshold impulse response, or a probe (Nikitin et al., 2003; Nikitin 
and Davidchack, 2003a,b). This threshold continuity of the output of a comparator is the 
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key to a truly analog representation of such a measurement. Examples of step and impulse 
responses of a continuous comparator are shown in figure 2. We further assume, for simplicity, 
that the probe is a unimodal even function, that is, /ad(D) has only a single maximum and 
ho(D) = h D (-D)< 

5 In practice, many different circuits can serve as comparators, since any continuous mono- 

tonic function with constant unequal horizontal asymptotes will produce the desired response 
under appropriate scaling and reflection. It may be simpler to implement a comparator de- 
scribed by an odd function T± D which relates to the response 

^ AD = A(2fkD-l), (D-l) 
10 where A is an arbitrary (nonzero) constant. For example, the voltage-current characteristic 
of a subthreshold transconductance amplifier (Mead, 1989; Urahama and Nagao, 1995) can 
be described by the hyperbolic tangent function, F AD = A tanh(D/AD) , and thus such an 
amplifier can serve as a continuous comparator. For specificity, this response function is used 
in the numerical examples of this disclosure. A practical implementation of the probe / AD 
15 corresponding to the comparator T± D can be conveniently accomplished as a finite difference 

Ud{d) w 5 -^6D- = TXW 1^° {D + 5D) ~ ^° {D - SD) 1 ' (D " 2) 

where 5D is a relatively small fraction of AD. 

Note that the terms Comparator' and 'discriminator' might be used synonymously in this 
disclosure. 

D-2 Analog rank filters (ARFs) 

20 

Consider the measuring process in which the difference between the threshold variable D and 
the scalar signal x(t) is passed through a comparator T± D , followed by a linear time averaging 
filter with a continuous impulse response w(t). The output of this system can be written as 

t) = w(t) * T^d [D - x(t)] , (D-3) 
25 where the asterisk denotes convolution. The physical interpretation of the function $(D,t) is 
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the (time dependent) cumulative distribution function of the signal x{t) in the moving time 
window w(t) (Nikitin and Davidchack, 2003b). In the limit of high resolution (small AD), 
equation (D-3) describes the 'ideal' distribution (Ferreira, 2001). Notice that $(D, t) is viewed 
as a function of two variables, threshold D and time £, and is continuous in both variables. 
5 The output of a quantile filter of order q in the moving time window w(t) is then given by 

the function D q (t) defined implicitly as 

*[D q {t),t]=q, 0<q<l. (D-4) 

Viewing the function $(£>, t) as a surface in the three-dimensional space (i, Z), $), we imme- 
diately have a geometric interpretation of D q (t) as that of a level (or contour) curve obtained 

10 from the intersection of the surface $ = ®(D,t) with the plane $ = q, as shown in figure 3. 
Based on this geometric interpretation, one can develop various explicit as well as feedback rep- 
resentations for analog rank filters, including such generalizations as L filters and a-trimmed 
mean filters (Nikitin and Davidchack, 2003b). 

Note that the terms analog 'rank', 'quantile', 'percentile', and 'order statistic' filters are 

15 often synonymous and might be used alternatively in this disclosure. 

D-2.1 Quantile filters for modulated densities 



Let us point out (see Nikitin and Davidchack, 2003a,b, for example) that various thresh- 
old densities can be viewed as different appearances of a general modulated threshold den- 
20 sity (MTD) 

a( n t\ ^)*{^)A.p-x(t)]} 

where K(t) is a unipolar modulating signal. Various choices of the modulating signal allow 
us to introduce different types of threshold densities and impose different conditions on these 
densities. For example, the simple amplitude density is given by the choice K(t) = const, 
25 and setting K(t) equal to \x(t)\ leads to the counting density. The significance of the defini- 
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tion of the time dependent counting (threshold crossing) density arises from the fact that it 
characterizes the rate of change in the analyzed signal, which is one of the most important 
characteristics of a dynamic system. Notice that the amplitude density is proportional to 
the time the signal spends in a vicinity of a certain threshold, while the counting density is 
5 proportional to the number of 'visits' to this vicinity by the signal. Figure 4 shows both the 
amplitude and counting densities computed for the fragment of a signal shown in the top panel 
of the figure. Note that the amplitude density has a sharp peak at every signal extremum, 
while the counting density has a much more regular shape. 

An expression for the quantile filter for a modulated density can be written as 

10 *£l w T (t)*{K(t)(q-^ D [D q -x(t)})} m fi x 

dt rw(t)*{K(t)f AD [D q -x(t))} ' KU ' D) 

and the physical interpretation of such a filter depends on the nature of the modulating signal. 

For example, a median filter in a rectangular moving window for K(t) = \x(t)\ f ADi [x(t)] yields 

Di/ 2 (t) such that half of the extrema of the signal x(t) in the window are below this threshold. 

Selected Acronyms and Where They First Appear 



15 AARF Adaptive Analog Rank Filter, page 27 

AARS Adaptive Analog Rank Selector, page 35 

ARTEMIS Adaptive Real-Time Embodiments for Multivariate Investigation of Signals, 

page 6 

AVATAR Analysis of Variables Through Analog Representation, page 4 

20 BASIS Bimodal Analog Sensor Interface System, page 40 

EARL Explicit Analog Rank Locator, page 37 

MTD Modulated Threshold Density, page 13 

SPART Single Point Analog Rank Tracker, page 33 
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Selected Notations and Where They First Appear 

9(x) Heaviside unit step function, page 7 

Fad, ?ad continuous comparator (discriminator), equation (D-l) on page 12 

D q {t), D q {t\T) output of quantile filter of order pages 13 and 41 

2?(D,x) threshold domain function, equation (1) on page 17 

\x\± positive/negative component of x, equation (7) on page 19 

5(x) Dirac delta function, equation (8) on page 19 

7£(D,£), 7l(t) instantaneous counting rate, page 19 and equation (8) on page 19 

i?(D,t), R(t) counting rate in moving window of time, page 19 and equation (17) on page 24 

delayed comparator, page 30 

!P££ averaging comparator, equation (27) on page 31 
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DETAILED DESCRIPTION OF THE INVENTION 

The Detailed Description of the Invention is organized as follows. 
Section 1 (p. 17) provides the definition of the threshold domain function and Threshold 
Domain Filtering, and explains its usage for feature extraction. 
5 Section 2 (p. 19) deals with quantification of crossings of threshold domain boundaries by 

means of Analog Counting. 

Section 3 (p. 22) introduces Multimodal Pulse Shaping as a way of embedding an incoming 
signal into a threshold space and thus enabling extraction of the features of interest by the 
Threshold Domain Filtering. Subsection 3.1 describes Analog Bimodal Coincidence (ABC) 
10 counting systems as an example of a real-time signal processing utilizing Threshold Domain 
Filtering in combination with Analog Counting and Multimodal Pulse Shaping. 

Section 4 (p. 26) presents various embodiments of Analog Rank Filters which can be used in 
ARTEMIS in order to reconcile the conflicting requirements of the robustness and adaptability 
of the control levels of the Threshold Domain Filtering. Subsection 4.1 describes the Adaptive 
15 Analog Rank Filters (AARFs) and Adaptive Analog Rank Selectors (AARSs), while §4.2 
introduces the Explicit Analog Rank Locators (EARLs). Subsection 4.3 describes the Bimodal 
Analog Sensor Interface System (BASIS) as an example of an analog signal processing module 
operatable as a combination of Threshold Domain Filtering, Analog Counting, and Analog 
Rank Filtering. 

20 As an additional illustration of ARTEMIS, §5 (p. 46) describes a technique and a circuit 
for generation of monoenergetic Poissonian pulse trains with adjustable rate and amplitude 
through a combination of Threshold Domain Filtering and Analog Counting. 
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1 Threshold domain filtering 

Threshold Domain Filtering is used for separation of the features of interest in a signal from the 
rest of the signal. In terms of a threshold domain, a 'feature of interest' is either a point inside of 
the domain, or a point on the boundary of the domain. In an electrical apparatus, e.g., a typical 

5 Threshold Domain Filter can be composed of (asynchronous) comparators and switches, where 
the comparators operate on the differences between the components of the incoming variable(s) 
and the corresponding components of the control variable(s). For example, for the domain 
defined as a product of two ideal comparators represented by the Heaviside unit step function 
V = 9[x(t) — D] 0[x(t)] (with two control levels, D and zero), a point inside (that 

10 is, V = 1) corresponds to a positive-slope signal of the magnitude greater than D, and 
the stationary points of x(t) above the threshold D can be associated with the points on the 
boundary of this domain. More generally, as used in the present invention, a Threshold Domain 
Filter is defined by its mathematical properties regardless of their physical implementation. 
Defining threshold domain Let us assume that a continuous signal y = y(a,i) depends 

15 on some spatial coordinates a and time t. Thus, in a vicinity of (a, t), this signal can be 
characterized by its value y(a, t) at this point along with its partial derivatives dy(a, t)/dai 
and dy(a,i)/dt at this point. These values (of the signal and its derivatives) can be viewed 
as coordinates of a point x = x(a, t) in a threshold space, where the vector x consists of the 
signal y and its various partial derivatives. A particular feature of interest can thus be defined 

20 as a certain region in the threshold space as follows. 

Let us describe an 'ideal' threshold domain by a (two-level) function £>(D,x) such that 

r 

X>(D, x) = qi if x is inside domain 

< * , (1) 

P(D, x) = q 2 otherwise 

where D is a vector of the control levels of the threshold filter. Without loss of generality, we 
can set qi = 1 and g 2 = 0. For example, in a physical space, an ideal cuboid with the edge 
25 lengths a, 6, and c, centered at (zo,yo 3 2o), can be represented by 
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s(w)=( [ 1 - 4 (£ ? )] ( [ 1 . 4 (» ? ) ! ] 8 [ 1 . 4 (i ? )'], 



(2) 



where we have assumed constant control levels and thus V is a function of x, y, and z only. 
Note that for a 'real', or 'fuzzy', domain the transition from q 2 to q\ happens monotonically 
over some finite interval (layer) of a characteristic thickness Aa. The transition to a 'real' 

5 threshold domain can be accomplished, for example, by replacing the ideal comparators given 
be the Heaviside step functions with the 'real' comparators, 6 — > T± D . 

Note that an arbitrary threshold domain can be represented by a combination (e.g., poly- 
nomial) of several threshold domains. For example, the cuboid given by equation (2) can be 
viewed as a product of six domains with plane boundaries, or as a product (intersection) of 

10 two domains given by the rectangular cylinders 



V xy (x } y) = e 

and 

V yz {y,z)=6 



£>(x, y, z) = V xy {x, y) V yz (y, z). 

15 Features of a signal In terms of a threshold domain, a 'feature 5 of a signal is either a point 
inside of the domain, or a point on the boundary of the domain. For example, for the domain 
V = 9[x(t) - D] 0[x(t)] (with two control levels, D and zero), a point inside (that is, V = 1) 
corresponds to a positive-slope signal of the magnitude greater than D, and a point on the 
boundary of this domain is a stationary point of x(t) above the threshold D. 

20 One should notice that only a small fraction of the signal's trajectory might fall inside of 
the threshold domain, and thus the duration of the feature might be only a small fraction of 
the total duration of the signal, especially if a feature is defined as a point on the domain's 
boundary. Therefore, it is impractical to continuously digitize the signal in order to extract 
the desired short-duration features. To resolve this, ARTEMIS utilizes an analog technique 

25 for extraction and quantification of the salient signal features. 
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2 Analog counting 



In its simplest form, analog counting consists of three steps: (1) time-differentiation, (2) rec- 
tification, and (3) integration. The result of step 2 (rectification) is the instantaneous count 
rate 7£(D,£), and step 3 (integration) outputs the count rate i?(D, t) in a moving window of 
time wit), i?(D, t) = w(t) * ft(D, t). 

Counting crossings of threshold domain boundaries The number of crossings of the bound- 
aries of a domain V by a point following the trajectory x(4) during the time interval [0, T] 
can be written as 



N-f 
Jo 

10 for the total number of crossings, or 



t 1 

N ± = 

Jo 



dt 



dt 



^2>P,x(t)] 



dt 



2>[D,x(t)] 



(5) 



(6) 



for the number of entries (+) or exits (— ). In equation (6), |x|± denotes positive/negative 
component of x, 



\x\ ± = \(\x\±x) 



(7) 



15 



Instantaneous count rates Note that the integrands in equations (5) and (6) represent the 
instantaneous rates of crossings of the domain boundaries, 



ft(r) = £<*(*-*<). 



(8) 



where 6(t) is the Dirac delta function, and U are the instances of the crossings. It should be 
20 easy to see that a number of other useful characteristics of the behavior of the signal inside 
the domain can be obtained based on the domain definition given by equation (1). 

Consider, for example, a threshold domain V in a physical space given by a product 
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(intersection) of two fields of view (e.g., solid angles) of two lidars 2 or cameras. When an 
object following the trajectory x(t) is in a field of view, the signal is T. Otherwise, it is '0'. 
Then the product of the signals from both lidars (cameras) is given by X>[D,x(i)], and the 
counting of the crossings of the domain boundaries by the object can be performed by an 
5 apparatus implementing equations (5) or (6). The following characteristics of the object's 
motion though the domain are also useful and easily obtained: The time spent inside the 
domain, 

t= / T dtP[D,x(*)], (9) 
Jo 

the distance traveled inside the domain, 

10 s= / r dt|x(t)|2>[D,x(t)], (10) 

Jo 

and the average speed inside the domain, 

fidt |x(t)|P[D,x(t)] 
t / 0 T diP[x(0] 



When using the 'real' comparators in a threshold filter and 'real' differentiators in an ana- 
15 log counter, the main property of the 'real' instantaneous rate TZ(t) is 

lim Hmft(t)=y>(t-ti), (12) 

t 

where AD and 5t are the width and the delay parameters of the comparators and differen- 
tiators, respectively. The property given by equation (12) determines the main uses of the 
instantaneous rate. For example, multiplication of the latter by a signal x(t) amounts to sam- 
20 pling this signal at the times of occurrence of the events U. Other temporal characteristics 
of the events can be constructed by time averaging various products of the signal with the 
instantaneous rate. 

Count rate in a moving window of time Count rate in a moving window of time w T (t) is 
obtained through the integration of the instantaneous rate by an integrator with an impulse 
25 response wr{t), namely as 

2 Here LEDAR is an acronym for "Light Detector And Ranger". 
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R(D, t) = wr(t) * ft(D, t) . (13) 

Although there is effectively no difference between averaging window functions which rise 
from zero to a peak and then fall again, boxcar averaging is deeply engraved in modern 

5 engineering, partially due to the ease of interpretation and numerical computations. Thus 
one of the requirements for counting with a non-boxcar window is that the results of such 
measurements are comparable with boxcar counting. As an example, let us consider averaging 
of the instantaneous rates by a sequence of n RC-integrators. For simplicity, let us assume 
that these integrators have identical time constants r = RC, and thus their combined impulse 

10 response is 

Wn{t) = J^T)^ e ~ t/Tm - (14) 

Comparability with a boxcar function of the width T can be achieved by equating the first 
two moments of the respective weighting functions. Thus a sequence of n RC-integrators with 
identical time constants r = ^T/s/3n will provide us with rate measurements corresponding 

15 to the time averaging with a rectangular moving window of width T. 3 

Figure 5 compares the moving window rates measured with the boxcar (thin solid line) 
and the 'triple-integrator' test function (n = 3 in equation (14), thick line) with r = T/6. The 
respective test functions are shown in the upper left corner of the figure. The gray band in 
the figure outlines the error interval in the rate measurements as the square root of the total 

20 number of counts in the time interval T per this interval. 

One of the obvious shortcomings of boxcar averaging is that it does not allow meaningful 
differentiation of counting rates, while knowledge of time derivatives of the event occurrence 
rate is important for all physical models where such rate is a time-dependent parameter. 
Indeed, the time derivative of the rate measured with a boxcar function of width T is simply 

25 T~ l times the difference between the 'original' instantaneous pulse train and this pulse train 

3 Of course one can design different criteria for equivalence of the boxcar weighting function and w(t). In 
our example we were simply looking for the width parameter of w(t) which would allow us to interpret the 
rate measurements with this function as 4 a number of events per time interval T\ 
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delayed by T, and such representation of the rate derivative hardly provides physical insights. 
On the other hand, the time derivative of the 'cascaded integrators' weighting function w n 
given by equation (14) is the bipolar pulse w n = r" 1 (u/ n -i — w n ), and thus the time derivative 
of the rate evaluated with w n (t) is a measure of the 'disbalance' of the rates within the 
5 moving window (positive for a 'front-loaded' sample, and negative otherwise). 

3 Multimodal pulse shaping 

In order to focus upon characteristics of interest, feature definition may require knowledge 
of the (partial) derivatives of the signal. For example, in order to count the extrema in a 

10 signal x(t), one needs to have access to the time derivative of the signal, x(t). A typical 
Multimodal Pulse Shaper in the present invention transforms at least one component of the 
incoming signal into at least two components such that one of these two components is a 
(partial) derivative of the other, and thus Multimodal Pulse Shaping can be used for embedding 
the incoming signal into a threshold space and enabling extraction of the features of interest 

15 by the Threshold Domain Filtering. 

Note, however, that differentiation performed by any physical differentiator is not accurate. 
For example, a time derivative of f(t) obtained by an RC differentiator is proportional to 
[e~ t/r 0(i)] */(£), where r = RC, not to f(t). Thus Multimodal Pulse Shaping does not attempt 
to straightforwardly differentiate the incoming signal. Instead, it processes an incoming signal 

20 in parallel channels to obtain the necessary relations between the components of the output 
signal. For example, if x(a, t) is a result of shaping the signal y(a, t) by the first channel of a 
pulse shaper with the impulse response /(a,i), 

x(a,i) = /(M)*y(a,t), (15) 
then the derivatives of x can be obtained as 

25 x(a, t) = / * y(a, t) and ^- x(a, t) = J£ * y(a, t) . (16) 
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Thus Multimodal Pulse Shaping will be achieved if the impulse responses of various channels 
in the pulse shaper relate as the respective derivatives of the impulse response of the first 
channel. 

Figure 6 shows an example of bimodal pulse shaping which can be used for both the 
5 amplitude and timing measurements of a short-duration event. An event of magnitude E { 
and arrival time U is passed through an RC pulse shaping network, producing a continuous 
signal x(t). The event can be fully characterized, e.g., by the first extremum of x(t), since 
the height of the extremum is proportional to and its position in time is delayed by a 
constant with respect to U. By replacing an RC integrator in the shaping network by an RC 
10 differentiator with the same time constant, one can obtain an accurate time derivative of x(t). 
Now the event can be associated with the inbound crossing of the boundary of the threshold 
domain V = 6[x(t) — D] Q[— x(t)] } where the threshold D is set at a positive value to eliminate 
the rest of the signal's stationary points. 

3.1 Example: Analog Bimodal Coincidence (ABC) counting systems 

15 Let us illustrate the usage of a threshold filter, in combination with multimodal pulse shap- 
ing and analog counting, in a signal processing module for a two-detector charged particle 
telescope. This module is an example of an Analog Bimodal Coincidence (ABC) counting 
system. 

In our approach, we relate the short-duration particle events to certain stationary points 
20 (e.g., local maxima) of a relatively slow analog signal. Those points can be accurately identified 
and characterized if the time derivative of the signal is available. Thus the essence of ABC 
counting systems is in the use of multiple signal characteristics - here a signal and its time 
derivative - and signals from multiple sensors in coincidence to achieve accuracy in both the 
amplitude and timing measurements while using low-speed, analog signal processing circuitry. 
25 This allows us to improve both the engineering aspects of the instrumentation and the quality 
of the scientific data. 

A simplified schematic of the module is shown in figure 7. A bimodal pulse shaping is used 
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to obtain an accurate time derivative of the signal from a detector. Comparators are used 
to obtain two-level signals with the transitions at appropriate threshold crossings (e.g., zero 
crossings for the derivative signal). Simple asynchronous analog switches are used to obtain 
the products of the comparators' outputs suitable for appropriate conditional and coincidence 

5 counting. The comparators and the analog switches constitute the threshold domain filter with 
the thresholds {^2}, and the grounds as the control levels. A-Counters are employed 

for counting the crossings of the threshold domains' boundaries. In its simplest form, an A- 
Counter is a differentiating circuit (such as a simple RC-differentiator) with a relatively small 
time constant (in order to keep the dead- time losses small), followed by a precision diode 

10 and an integrator with a large time constant (at least an order of magnitude larger than the 
inverted smallest rate to be measured). A TOF selector employs an additional pulse shaping 
amplifier, and a pair of comparators with the levels corresponding to the smallest and the 
largest time of flight. 

Bimodal pulse shaping and instantaneous rate of signal's maxima When the time derivative 
15 of a signal is available, we can relate the particle events to local maxima of the signal and 
accurately identify these events. Thus bimodal pulse shaping is the key to the high timing 
accuracy of the module. As shown in figure 7, a bimodal pulse shaping unit outputs two 
signals, where the second signal is proportional to an accurate time derivative of the first 
output. The rate R(t), in the moving window of time wr(t), of a signal's maxima above the 
20 threshold D can be expressed as 



R(t) = wr(t) * 



±{0[z{t)-D]9[-±m 



(17) 



+ 

where denotes the positive part of y (see equation (7)), 0 is the Heaviside unit step 
function, and the asterisk denotes convolution. Equation (17) represents an idealization of 
the measuring scheme consisting of the following steps: (i) the first output of the bimodal 
25 pulse shaping unit is passed through a comparator set at level D, and the second output - 
through a comparator set at zero level; (ii) the product of the outputs of the comparators is 
differentiated, (iii) rectified by a (precision) diode, and (iv) integrated on a time scale T (by 
an integrator with the impulse response wr{t)). Note that steps (ii) through (iv) represent 
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passing the product of the comparators through an A-Counter. Also note that the output of 
step (iii) is the instantaneous rate of the signal's maxima above threshold D. 
Basic coincidence counting For basic coincidence counting, the coincident rate R c {t) can be 
written as 

" d 



R c (t) = w T (t) * 



df {0[xx(t) - A] *[-±i(t)l 0[x 2 (t) - D 2 ]} 



(18) 



where the notations are as in equation (17). One can see that equation (18) differs from 
equation (17) only by an additional term in the product of the comparators' outputs. 
Transition to realistic model of measurements It can be easily seen that equations (17) 
and (18) do not correctly represent any practical measuring scheme implementable in hard- 

10 ware. For example, both equations contain derivatives of discontinuous Heaviside functions, 
and thus instantaneous rates axe expressed through singular Dirac J-functions. To make a 
transition from an ideal measurement scheme to a more realistic model, we replace the Heavi- 
side step functions by 'real' discriminators (6(x) — ► ot 6t {t)*T AD (x), where a St (t) is a continuous 
kernel such that J^dt a st (t) = 1), and perform differentiation through a continuous kernel 

15 (di • • - ~* * • • ■)> etc - We choose appropriate functional representations of T± D , aj 5t (i), 
etc., for various elements of a schematic, and also add appropriate noise sources such as 
thermal noise at all intermediate measuring steps. Figures 8 and 9 illustrate such realistic 
measurements of instantaneous rates of extrema and coincident maxima, respectively. Notice 
that, in both figures, an event is represented by a narrow peak of a prespecified area in the 

20 instantaneous rates. 

Time-of-flight (TOF) constrained measurements The time-of-flight constrained coincident 
rate can be expressed, for times of flight larger than At, as 



1Z C = w T (t) * 



^0[h*\V 12 -V 21 \ + -Z At ] 



(19) 



+ 

where h is some (unipolar or bipolar) impulse response function, Z& t is a threshold level 
25 corresponding to the TOF equal to At, and D {j = \0[xi(t) - A] 0 [-£<(£)] 6 [x/yt) - Dj]\ + . 
Thus a TOF selector (see figure 7) will consist of a pulse shaping amplifier with an impulse 
response h, and a differential comparator. Figure 10(a) illustrates coincident counting 
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according to equation (19), and figure 10(b) provides an example of using a realistic model 
of the TOF measurements, with functional representations of the elements of the schematic 
corresponding to commercially-available, off-the-shelf (COTS) components. As can be seen 
in the figure, the performance of the system is not significantly degraded by the transition 
5 from an idealized to a more realistic model. 

4 Analog Rank Filtering 

In ARTEMIS, Analog Rank Filtering can be used for establishing and maintaining the analog 
control levels of the Threshold Domain Filtering. It ensures the adaptivity of the Threshold 

10 Domain Filtering to changes in the measurement conditions (e.g., due to nonstationarity of the 
signal or instrument drift), and thus the optimal separation of the features of interest from the 
rest of the signal. For example, the threshold level D in the domain V = 0[x(t) — D] 0[x(t)] can 
be established by means of Analog Rank Filtering to separate the stationary points of interest 
from those caused by noise. Note that the Analog Rank Filtering outputs the control levels 

15 indicative of the salient properties of the input signal(s), and thus can be used as a stand-alone 
embodiment of ARTEMIS for adaptive real-time signal conditioning, processing, analysis, 
quantification, comparison, and control, and for detection, quantification, and prediction of 
changes in signals. 

Creating and maintaining baseline and analog control levels by analog rank filters Analog 
20 rank filters can be used to establish various control levels (reference thresholds) for the thresh- 
old filter. When used in ARTEMIS, rank-based filters allows us to reconcile, based on the 
rank filters' insensitivity to outliers, the conflicting requirements of the robustness and adapt- 
ability of the control levels of the Threshold Domain Filtering. In addition, the control levels 
created by Analog Rank Filters axe themselves indicative of the salient properties of the input 
25 signai(s). 
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4.1 Adaptive Analog Rank Filters (AARFs) 

Rank filter in RC window When the time averaging filter in equation (D-3) is an RC integra- 
tor (RC = r), the differential equation for the output D q (t) of a rank filter takes an especially 
simple form and can be written as 

dD q __ A(2q-1)-F AD [D q (t)-x(t]\ 



(20) 



s=t 



dt 2ATh r {s)*f AD [D q (t)-x(s)} 

where h T (t) = 6{t) exp [— \ — lnr). 4 The solution of this equation is ensured to rapidly 
converge to D q (t) of the chosen quantile order q regardless of the initial condition (Nikitin 
and Davidchack, 2003b). Note also that the continuity of the comparator is essential for the 
right-hand side of equation (20) to be well behaved. 

10 The main obstacle to a straightforward analog implementation of the filter given by equa- 
tion (20) is that the convolution integral in the denominator of the right-hand side needs to be 
re-evaluated (updated) for each new value of D q . If we wish to implement an analog rank filter 
in a simple feedback circuit, then we should replace the right-hand side of equation (20) by 
an approximation which can be easily evaluated by such a circuit. Of course, one can employ 

15 a great variety of such approximations (Bleistein and Handelsman, 1986, for example), whose 
suitability will depend on a particular goal. A very simple approximation becomes available 
in the limit of sufficiently small r, since then we can replace h T (s) * f AD [D q (t) — x(s)] | by 
h T (t) * /az> [D q (t) — x(t)] in equation (20). As was shown by Nikitin and Davidchack (2003b), 
this simple approximation can still be used for an arbitrary time window w(t), if we repre- 

20 sent w(t) as a weighted sum of many RC integrators with small r. However, this approximation 
fails when the threshold resolution is small (e.g., when AD < \h r (t) * x(t)\ r), and thus cannot 
be used in real-time processing of non-stationary signals. 



4 In more explicit notation, the convolution integral in the denominator of equation (20) can be written as 
h T (s) * /ad [D q (t) - *(*)] \ s=t = \f_te exp {^rj Ud [D q {t) ~ *«] ■ 
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Adaptive approximation of a feedback rank filter in an arbitrary time window A rank filter 
in a boxcar moving time window B T (t) = [8(t) — 9(t - T)] /T is of a particular interest, since 
it is the most commonly used window in digital rank filters. The output D q of an analog 
rank filter in this window is implicitly defined as Br{t) * T± D [D q — x(t)] = q. To construct 
5 an approximation for this filter suitable for implementation in an analog feedback circuit, we 
first approximate the boxcar window B T {t) by the following moving window wpf(t): b 

%(t)4EMt-2fcr), (21) 
iV k=0 

where r = T/(2N). The first moments of the weighting functions w N (t) and B T (t) are 

identical, and the ratio of their respective second moments is yjl + 2/N 2 » 1+1 /N 2 . The other 

10 moments of the time window wjv(0 also converge rapidly, as N increases, to the respective 

moments of B T {t), which justifies the approximation of equation (21). 

Now, the output of a rank filter in such a window can be approximated as discussed earlier, 

namely as (Nikitin and Davidchack, 2003b) 6 

dD, _ A N(2g - 1) - E^o 1 [D q (t) - x(t - 2kr)) 
dt ~ 2Arh T (t)*Z^ h D [D q (t)-x(t-2kr)] J 1 ; 

15 where r = T/(2N). Note that the accuracy of this approximation is contingent on the re- 
quirement that AZ? > \h T (t) * x(t)\ r. This means that, if we wish to have a simple analog 
circuit and keep N relatively small, we must choose AD sufficiently large for the approxima- 
tion to remain accurate. On the other hand, we would like to maintain high resolution of the 
acquisition system, that is, to keep AD small. 

20 In order to reconcile these conflicting requirements, we propose to use an adaptive approxi- 
mation, which reduces the resolution only when necessary. This can be achieved, for example, 
by using equation (D-2) and rewriting the threshold derivative of h T (t) * T± D [D q — x(t)] as 



5 Since a moving time window is always a part of a convolution integral, the approximation is understood 

in the sense that Br(t) * g(t) « w^it) * g{t), where g(t) is a smooth function. 

6 An explicit expression for the convolution integral hr(t) * f AD [D q (t) — x(t — 2kr)] is 



h T (t) * /ad \D q {t) - x(t - 2kr)] = i J ds exp (j^r^j /ad [D q (s) - x(s - 



2kr)\ 
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Mg * {^g [D q+ - x(t)} - p q _ - 
M«) * ho [D q - x(t)\ * L 4A(D 9+ -!>,-) ^ ' (23) 

where D q ± is the output of a rank filter of the quantile order q ± 5q, 8q q. In essence, the 
approximation of equation (23) amounts to decreasing the resolution of the acquisition system 
only when the amplitude distribution of the signal broadens, while otherwise retaining high 
resolution. 

Combining equations (21-23), we arrive at the following representation of an adaptive 
approximation to a feedback rank filter in a boxcar time window of width T: 

D g (t) = \ [D q+ (t) + D q .{t)\ 

n , A A N(2q - 1 + 25q) - Egg jkg [D g+ (t) ~ x(t - 2kr)) 8D q (t) 

D < +it) = h T {t)*6TUt) ~ • ( 24 ) 

^ ( ) A N(2q - 1 - 25q) - E^'q 1 [D q - (t) - x(t - 2fcr)] 5D q (t) 
" K{t) * 5f AD (t) t 

where 5D q (t) = D q+ (t) - D q .(t) and 

SfUt) = E l D Mt) ~ <t ~ 2*r)] - JF AD [£>,_(«) - X (t - 2kr)}} . (25) 

fc=0 

This approximation preserves its validity for high resolution comparators (small AD), and its 
output converges, as N increases, to the output of the 'exact' rank filter in the boxcar time 
window B T {t). Unlike the currently known approaches (see, for example, Urahama and Nagao 
1995; Opris 1996), the analog rank filters enabled through equation (24) are not constrained by 

15 linear convergence and allow real time implementation on an arbitrary timescale, thus enabling 
high speed real time rank filtering by analog means. The accuracy of this approximation is 
best described in terms of the error in the quantile q. That is, the output D q (t) can be viewed 
as bounded by the outputs of the 'exact' rank filter for different quantiles q ± Aq. When AD 
and 8q in equation (24) are small, the error range Aq is of order 1/N. 

20 Note that, even though equation (24) represents a feedback implementation of a rank filter, 
it is stable with respect to the quantile values q. In other words, the solution of this equation 
will rapidly converge to the 'true' value of D q (t) regardless of the initial condition, and the 
time of convergence within the resolution of the filter AD for any initial condition will be just a 
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small fraction of r. This convergence property is what makes the implementation represented 
by equation (24) suitable for a real time operation on an arbitrary timescale. 
Implementation of AARFs in analog feedback circuits Figure 11 illustrates implementation 
of an adaptive real time rank filter given by equation (24) in an analog feedback circuit. One 
5 skilled in the art will recognize that this circuit is a simplified embodiment of a more general 
AARF depicted in figure 12. 

Generalized description of AARFs As shown in figure 12, an input variable x(t) and a plu- 
rality of feedbacks of Offset Rank Filtered Variables {D qi (t)} are passed through a plural- 
ity of (delayed) comparators forming a plurality of outputs of the comparators {^^(i)} = 

10 {Fto [D qi (t),x(t)]}. (Please note that in this and further figures a double line in a diagram 
indicates a plurality of signals.) Said plurality of the outputs of the comparators {T* eX } is 
used to form (i) a plurality {A(2qi — 1) — i^ e1 } of differences between said outputs of the 
comparators and the respective Offset Quantile Parameters of said Offset Rank Filtered Vari- 
ables, and (ii) a weighted difference 5F AD (t) = T,i<XiF? e \ where Ei^i = 0, of said outputs 

15 of the comparators. Said weighted difference 5F AD (t) of the outputs of the comparators is 
passed through a time averaging amplifier, forming a Density Function h T (t) * 6F AD (t). The 
plurality of the feedbacks of the Offset Rank Filtered Variables {D qi (t)} is used to form a 
weighted difference 8D q (t) — J2iPiD qi (t), where 5^ fk = 0, of said feedbacks. Each difference 
A(2qi — 1) — Tf^ between the outputs of the comparators and the respective Offset Quan- 

20 tile Parameters of the Offset Rank Filtered Variables is multiplied by a ratio of the weighted 
difference SD q (t) of the feedbacks of the Offset Rank Filtered Variables and the Density Func- 
tion h r (t) * 8T AD (t), forming a plurality of time derivatives of Offset Rank Filtered Variables 
{D qi (t)}. Said plurality of the time derivatives {D qi {t)} is integrated to produce the plurality 
of the Offset Rank Filtered Variables {D qi (t)}> The plurality of the Offset Rank Filtered 

25 Variables {D qi (t)} is then used to form an output Rank Filtered Variable D q (t) as a weighted 
average £ f WiD Qi (t), ttfi = 1, of said Offset Rank Filtered Variables. 

As an example, figure 13 provides a simplified diagram of a 3-comparator implementation 
of AARF. In this example, the offset quantile orders are q<i = q, q x = q - dq, and <? 3 = q + dq, 
and the weights for weighted average and differences are: w 2 = 1, w x = w 3 — 0, a 2 = fa = 0, 
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and c*i = Pi = -a 3 = ~P$ — -1. 

Note that both the input and output of an AARF are continuous signals. The width of the 

moving window and the quantile order are continuous parameters as well, and such continuity 

can be utilized in various analog control systems. The adaptivity of the approximation allows 
5 us to maintain a high resolution of the comparators regardless of the properties of the input 

signal, which enables the usage of this filter for nonstationary signals. 

Also, let us point out that the equations describing this filter are also suitable for numerical 

computations, especially when the number of data points within the moving window is large. 

A simple forward Euler method is fully adequate for integrating these equations, and the 
10 numerical convolution with an RC impulse response function requires remembering only one 

previous value. Thus numerical algorithms based on these equations have the advantages of 

both high speed and low memory requirements. 

Delayed comparators In our description of AARFs we have assumed that the comparators 
are the delayed comparators with the outputs represented by the moving averages 

15 Tit [D q (t),x(t)) = £ W) ~ *(* - A <*)1 , (26) 

where it;* are positive weights such that Y,k w k = 1, and it can be assumed, without loss of 
generality, that Ato = 0. Obviously, when N = 1, a delayed comparator is just a simple 
two-level comparator. Figure 14 illustrates a principle schematic of a delayed comparator. 
Averaging comparators In the description of Adaptive Analog Rank Selectors further in this 
20 disclosure we will use another type of a comparator, which we refer to as an averaging com- 
parator. Unlike a delayed comparator which takes a threshold level and a scalar signal as 
inputs, the inputs of an averaging comparator are a threshold level D and a plurality of input 
signals {xi(t)}, i = 1, ... ,7V. The output of an averaging comparator is then given by the 
expression 

25 J*S [DM {*#)}] = £>, f. D \D q (t) - x*(t)] , (27) 

1=1 

where Wi are positive weights such that Wi = 1. Figure 15 illustrates a principle schematic 
of an averaging comparator. 
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Figure 16 compares the performance of the analog rank filter given by equation (24) to 
that of the 'exact' quantile filter in a boxcar moving window of width T. In this example, the 
quantile interval Sq is chosen as Sq = 10~ 2 (1%). The continuous input signal x(t) (shown 
by the solid dark gray line) is emulated as a high resolution time series (2 x 10 3 points per 
interval T). The 'exact' outputs of a boxcar window rank filter are shown by the dashed lines, 
and their deviations within the ±Aq intervals are shown by the gray bands. The respective 
outputs of the approximation given by equation (24) are shown by the solid black lines. The 
width parameter AD of the comparators, the width T of the boxcar time window, the quantile 
order q, and the number TV of exponential kernels in the approximation are indicated in the 
figure. 

The (instantaneous) accuracy of the approximation given by equation (24) decreases when 
the input signal x(t) undergoes a large (in terms of the resolution parameter AD) monotonic 
change over a time interval of order r. The main effect of such a 'sudden jump' in the input 
signal is to delay the output D q (t) relative to the output of the respective 'exact' filter. This 
delay is shown as At in the lower left portion of the upper panel, where the input signal is a 
square pulse. This timing error At is inversely proportional to the number N of the kernels 
in the approximation. The accuracy of the approximation can also be described in terms of 
the amplitude error. As can be seen in figure 16, the residual oscillations of the outputs of the 
analog filter occur within the q ± l/(2N) interval around the respective outputs of the 'exact' 
filter (that is, within the width of the gray bands in the figure). 
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Establishing internal reference signal (baseline and analog control levels) As stated earlier, a 
primary use of Analog Rank Filtering in ARTEMIS is establishing and maintaining the analog 
control levels of the Threshold Domain Filtering, which ensures the adaptivity of the Threshold 
Domain Filtering to changes in the measurement conditions, and thus the optimal separation 
5 of the features of interest from the rest of the signal. Such robust control levels can be 
established, for example, by filtering the components of the signal with a Linear Combination 
of Analog Order Statistics Filters operable on a given timescale. 

Example: "Trimean' reference. Figure 17 provides an example of using an internal reference 
(baseline) for separating signal from noise, and illustrates a technique for establishing a refer- 

10 ence baseline as a linear combination of quartile outputs (i.e., q = 1/4, q = 1/2, and q = 3/4) 
of AARFs. In this example, the features of interest are tall pulses protruding from a noisy 
background. For example, one would want to count the number of such pulses, while ignoring 
the smaller pulses due to noise. This can be accomplished by choosing a reference baseline 
such that most of the pulses of interest peak above this baseline, while the accidental crossings 

15 of the baseline by noise are rare. A good choice for a baseline thus would be a moving average 
of the noise plus several standard deviations of the noise in the same moving window of time 
(a 'variance 5 baseline, gray lines in the figure). However, the presence of the high-amplitude 
pulses of the 'useful' signal will significantly disturb such a baseline. Instead, one can create a 
baseline by using a linear combination of the outputs of AARFs for different quantile orders 

20 (e.g., for the quartiles q = 1/4, q — 1/2, and q = 3/4 - 'quartile' baseline, dashed lines in 
the figure). As shown in the upper panel of the figure, in the absence of the signal of interest 
the baselines created by both techniques are essentially equivalent. However, as shown in the 
lower panel of the figure, in the presence of tall pulses the 'variance' baseline is significantly 
disturbed and fails to separate the noise from the signal, while the 'quartile' baseline remains 

25 virtually unaffected by the addition of these pulses. In both panels, the distance between the 
time ticks is equal to the width of the moving time window. 

Figure IS illustrates a principle diagram of a circuit for establishing a baseline as a linear 
combination of the quartile outputs (q = 1/4, q = 1/2, and q = 3/4) of AARFs. One skilled 
in the art will recognize that a variety of other linear combinations of outputs of AARFs of 
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different quantile orders can be used for establishing and maintaining the analog control levels 
of the Threshold Domain Filtering. 

Single Point Analog Rank Tracker (SPART) The approximation of equation (24) preserves its 
validity for high resolution comparators (small AD), and its output converges, as N increases, 

5 to the output of the 'exact 5 rank filter in the boxcar time window B T (t). However, even a single- 
point approximation (N = 1 in equation (24), i.e., simple rather than delayed comparators in 
AARF) can be fully adequate for creating and maintaining the baseline and analog control 
levels in analog counting systems, since such a simplified implementation preserves the essential 
properties of the 'exact 5 rank filter needed for this purpose. We shall call this version of an 

10 AARF the 'Single Point Analog Rank Tracker 5 , or SPART. 
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10 



Adaptive Analog Rank Selectors (AARSs) While an AARF operates on a single scalar input 
signal x(t) and outputs a qth quantile D g (t) of the input signal in a moving window of time, 
an AARS operates on a plurality of input signals {^(i)}, i — 1, . . . , N, and outputs ('selects') 
an instantaneous qth quantile D q (t) (in general, a weighted quantile) of the plurality of the 
input signals. Such transition from an AARF to an AARS can be achieved by replacing the 
delayed comparators in an AARF by averaging comparators. For example, a 2-comparator 
AARS can be represented by the following equation: 

D q (t) = \ [D q+ (t) + D q .(t)\ 

n m AN{2q - 1 + 25q) - JF AP [D q+ (t) - XjQ] SD q (t) 

D < +{t) = M*)*«^0 ~ ' (28) 

• = AN(2q - 1 - 26q) - Eli [£>,-(*) ~ Xj(t)] 5D q (t) 
" K ' h T {t) * 6^ D (t) r 

where 8D q (t) = D q+ (t) - D q _(t) and 

5TUt) = E {F*d [D q+ (t) ~ Xi(t)) - f ifl [D q .(t) - Xi (t)]} . (29) 

i=l 

Note that the time of convergence (or time of rank selection) is proportional to the time 
constant r = RC of the RC integrator, and thus can be made sufficiently small for a true real 
time operation of an AARS. Figure 19 illustrates a principle schematic of an Adaptive Analog 
Rank Selector given by equation (28) in an analog feedback circuit. One skilled in the art 
15 will recognize that this circuit is a simplified embodiment of a more general AARS depicted 
in figure 20. 

Generalized description of AARSs As shown in figure 20, a plurality input variables {xj(i)}, 
j = 1 ? . . . ? iV, and a plurality of feedbacks of Offset Rank Selected Variables {D qi (t)} are 
passed through a plurality of averaging comparators forming a plurality of outputs of the 
20 comparators {^? ve (t)} = {f™ [D qi (t), x^t)]} = {E^Li Vj ?ad [D qi (t) - Xj(t)]}. Said plurality 
of the outputs of the comparators {Jf ve } is used to form (i) a plurality {^4(2^ - 1) - Tf ve } 
of differences between said outputs of the comparators and the respective Offset Quantile 
Parameters of said Offset Rank Selected Variables, and (ii) a weighted difference 5F AD (t) = 
X)i^t^? ve j where = 0, of said outputs of the comparators. Said weighted difference 
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S^Ao(t) of the outputs of the comparators is passed through a time averaging amplifier, 
forming a Density Function h T (t) * 8T± D (t). The plurality of the feedbacks of the Offset 
Rank Selected Variables {D q .(t)} is used to form a weighted difference 5D q (t) = YliPiD qi (t), 
where # = 0, of said feedbacks. Each difference A(2qi — 1) — ^ ve between the outputs of 

5 the comparators and the respective Offset Quantile Parameters of the Offset Rank Selected 
Variables is multiplied by a ratio of the weighted difference 8D q (t) of the feedbacks of the Offset 
Rank Selected Variables and the Density Function hr(t)*5F AD (t) } forming a plurality of time 
derivatives of Offset Rank Selected Variables {D Qi (t)}. Said plurality of the time derivatives 
{D qi (t)} is integrated to produce the plurality of the Offset Rank Selected Variables {D q .(t)}. 

10 The plurality of the Offset Rank Selected Variables {D qi (t)} is then used to form an output 
Rank Selected Variable D q (t) as a weighted average w*jD 9i (£)> DiWt = 1, of said Offset 
Rank Selected Variables. 

Adaptive Analog Rank Selectors are well suited for analysis and conditioning of spatially- 
extended objects such as multidimensional images. For example, a plurality of input signals 

15 can be the plurality of the signals from a vicinity around the spatial point of interest, 
and the weights {vj} can correspond to the weights of a spatial averaging kernel. This 
enables us to design highly efficient real-time analog rank filters for removing dynamic as 
well as static impulse noise from an image, as illustrated in figure 21 for a two-dimensional 
monochrome image. In this example, a median filter (q = 1/2) according to equation (28) is 

20 used. 7 Panel (a) shows the original (uncorrupted) image. Panel (b) shows the snapshots, at 
different times, of the noisy image and the respective outputs of the filter. In this example, 
approximately 4/5 of the pixels of the original image are affected by a bipolar non-Gaussian 
random noise at any given time. Panel (c) provides an example of removing the static noise 
(1/3 of the pixels of the original image are affected). This example also illustrates the fact 

25 that the characteristic time of convergence of the filter based on equation (28) is only a small 

fraction of the time constant r = RC of the RC integrator, which makes this circuit suitable 

for a truly real-time operation. This fast convergence is a consequence of the fact that the 

7 In general, the quantile order of the filter should be chosen as q = $ n (0), where $ n is the amplitude 
distribution of the noise (either measured or known a priori). In the example of this section, $ n (0) = 1/2. 
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speed of convergence is inversely proportional to the density function h T (t) * 8T± D (t). 

4.2 Explicit Analog Rank Locators (EARLs) 

Explicit expression for an analog quantile filter Note that a differential equation is not the 
5 only possible embodiment of an analog quantile filter. Other means of locating the level lines 
of the threshold distribution function can be developed based on the geometric interpretation 
discussed in §D-2. For example, one can start by using the sifting property of the Dirac 8- 
function to write D q (t) as 

D q (t) = f°°dD D S[D - D q {t)] (30) 

J —CO 

10 for all t. Then, recalling that D q {t) is a root of the function <f>(D,t) — q and that, by con- 
struction, there is only one such root for any given time t, we can replace the 5-function of 
thresholds with that of the distribution function values as follows: 

D q (t) = f°°dD D <f>(D, t) S [*(£>, t) - q) . (31) 

J -co 

Here we have used the following property of the Dirac 5-function (see Davydov, 1988, p. 610, 
15 eq. (A 15), for example): 

%-/WI=E^f. (32) 

where \f(xi)\ is the absolute value of the derivative of f(x) at £», and the sum goes over all 
Xi such that /(£*) = a. We have also used the fact that <f>{D, t) > 0. 

The final step in deriving a practically useful realization of the quantile filter is to replace 
20 the 5-function of the ideal measurement process with a finite- width pulse function g Aq of the 
real measurement process, namely 

D q (t) = [°°dD D cf>(D, t) g* q t) - q] , (33) 

J —CO 

where Aq is the characteristic width of the pulse. That is, we replace the 5-function with a 
continuous function of finite width and height. This replacement is justified by the observation 
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made earlier: it is impossible to construct a physical device with an impulse response expressed 
by the ^-function, and thus an adequate description of any real measurement must use the 
actual response function of the acquisition system instead of the ^-function approximation. 
We shall call an analog rank filter given by equation (33) the Explicit Analog Rank Locator 
5 (EARL). 

Analog L filters and q-trimmed mean filters It is worth pointing out the generalization of 
analog quantile filters which follows from equation (31). In the context of digital filters, this 
generalization corresponds to the L filters described by Bovik et al. (1983). 

Indeed, we can write a linear combination of the outputs of various quantile filters as 

10 D L {t) = Jdq W L (q) D q (t) = fdq W L (q) JjLD D <f>(D, t) S [$(D, t) - q] 

= f°dD D t) W L [*(Z>, t)] , (34) 

where Wl is some (normalized) weighting function. Note that the difference between equa- 
tions (34) and (33) is in replacing the narrow pulse function g± q in (33) by an arbitrary 
weighting function Wl. 

A particular choice of W L in (34) as the rectangular (boxcar) probe of width 1 - 2a, 
15 centered at 1/2, will correspond to the digital a-trimmed mean filters described by Bendat 
(1998): 

TOO 

D a {t) = dDD 4>{D, t) b Q [$(£>, t)] , 0 < a < 1/2 , (35) 

J — OO 

where b a (x) = [0(x - a) - 9(x — 1 + a)]. When a — 0, equation (35) describes the 
running mean filter, D a= o(t) = x(t), and in the limit a — > 1/2 it describes the median filter, 

20 \im a ^ l/2 D a (t) ^ D m {t). 

Dealing with improper integration: Adaptive EARL The main practical shortcoming of the 
filter given by equation (33) is the improper integral with respect to threshold. This difficulty, 
however, can be overcome by a variety of ways. 

For example, we can use the fact that rank is not affected by a monotonic trans- 

25 formation. That is, if D q is the gth quantile of the distribution w T (t) * 0[D — x(t)] 
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(that is, w r (t) * 8[D q — x(t)] = q), then F{D q ) is the gth quantile of the distribution 
w T {t)*0{F{D)-F[x{t)\Y 

w r (t) * 9 {HD q ) - T [x(t)]} = q , (36) 
where J 7 ^) is a monotonically increasing function of f . 

Now let us choose £ = J 7 ^) as the response of a real comparator, ^(0 = ^* 2 (£ ~ 
where /xi is indicative of the mean value of x(t) in a moving window w T of the width T much 
greater than r, and the width parameter fj, 2 is indicative of the signal's deviation around fa 
(on a similar time interval). For example, 

MiW =**(*)**(*) • ( 3? ) 

/*.(*) = V 2 W*)*^(*)-M?W1 
10 Then an equation for the adaptive explicit analog rank locator can be rewritten as 

D q (t) = fa(t) + T~l yW X 4>(X, t) g± q [*(*, *)"«]} i (38) 



where 



and 



^ ( *' * } = w r (t)*K(t) ' (39) 



**<>- ^-™**-™ . (40) 

Note that the improper integral of equation (33) has become an integral over the finite interval 
[0, 1], where the variable of integration is a dimensionless variable x- 

Figure 22 illustrates the performance of adaptive EARLs operating as amplitude (panel b)) 
and counting (panel (c)) rank filters in comparison with the 'exact' outputs of the respective 
20 analog rank filters given by equation (D-6). 

Discrete-Threshold Approximation to Adaptive EARL Given a monotonic array of threshold 
values between zero and unity, the integral in equation (38) can be evaluated in finite differ- 
ences leading to a discrete-threshold approximation to adaptive EARL as follows: 
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D q (t) = to(t) + f£(D q ), (41) 

where ~D q is the root of <I(Z), t) = q. For example, D q — \{D jx + D j2 ), where 0 < A < 1 is a 
monotonic array of threshold values, A < A+i> an d ji and j 2 are such that ^(Dj^t) < q < 
^(Dj 1+U t) and $(D j2 ,t) < q < $(Dj 2 + x ,t). Note that a binary search, as well as more effec- 
5 tive methods, can be used for the root finding, and thus the discrete-threshold approximation 
to adaptive EARL can have significant advantages over the state-of-art numerical algorithms 
for rank filtering, especially when operating on large time scales. 

Discrete-Threshold Approximation to AARF It is worth pointing out that the invariance 
of rank to a monotonic transformation allows us to define the following discrete-threshold 
10 approximation to an adaptive analog rank filter: 

t = K(t) {g-^[D q (t)-m}} 

where D k (t) = SDk(t) = 5Dmnt(D q (t)/6D)? and 8D < AD. 

4.3 Example: Bimodal Analog Sensor Interface System (BASIS) 



15 BASIS constitutes an analog signal processing module, initially intended to be coupled with 
a photon counting sensor such as a photomultiplier tube (PMT). The resulting integrated 
photodetection unit allows fast and sensitive measurements in a wide range of light intensi- 
ties, with adaptive automatic transition from counting individual photons to the continuous 
(current) mode of operation. When a BASIS circuit is used as an external signal processing 

20 unit of a photosensor, its output Rout(t) is a continuous signal for both photon counting and 
current modes, with a magnitude proportional to the rate of incident photons. This signal can 
be, for example, used directly in analog or digital measuring and/or control systems, differ- 
entiated (thus producing continuous time derivative of the incident photon rate), or digitally 
sampled for subsequent transmission and/or storage. Thus, BASIS converts the raw output of 
8 The nearest integer function nint(x) is defined as the integer part of x + ^, nint(rr) = [x + |J. 
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a photosensor to a form suitable for use in continuous action light and radiation measurements. 

The functionality of the BASIS is enabled through the integration of three main com- 
ponents: (1) Analog Counting Systems (ACS), (2) Adaptive Analog Rank Filters (AARF), 
and (3) Saturation Rate Monitors (SRM), as described further. The BASIS system pro- 

5 vides several significant advantages with respect to the current state-of-art signal processing 
of photosignals. Probably the most important advantage is that, by seamlessly merging the 
counting and current mode regimes of a photosensor, the output of the BASIS system has a 
contiguous dynamic range extended by 20-30 dB. This technical enhancement translates into 
important commercial advantages. For example, the extension of the maximum rate of the 

10 photon counting mode of a PMT by 20 dB can be used for a tenfold increase in sensitivity or 
speed of detection. Since sensitivity and speed of light detecting units is often the bottleneck 
of many instruments, this increase will result in upgrading the class of equipment at a fraction 
of the normal cost of such an upgrade. 

In addition, the analog implementation of the current mode regime reduces the overall 

15 power consumption of the detector. These capabilities will benefit applications dealing with 
light intensities significantly changing in time, and where autonomous low-power operation is 
a must. One particular example of such an application is a high sensitivity handheld radiation 
detection system that could be powered with a small battery. Such a compact detector could 
be used by United States customs agents to search for nuclear materials entering the country. 

20 Principal components (modules) of BASIS As shown in figure 23, the principal components 
(modules) of the BASIS can be identified as (I) Rank Filtering (or Baseline) Module, (II) Ana- 
log Counting Module (ACM), (III) the Saturated Rate Monitor (SRM), and (IV) Integrated 
Output Module. A brief description of these modules is as follows. 

Rank Filtering (or Baseline) Module As shown in figure 23, the Baseline Module outputs the 
25 rank-filtered signal D q (t]T), which is the qth quantile of the signal x(t) in a moving time 
window of characteristic width T. The rank filtering is accomplished by means of an Adaptive 
Analog Rank filter (AARF) (see §4.1), or its single-point version referred to as a Single Point 
Analog Rank Tracker (SPART) (ibidem). AARFs, due to their insensitivity to outliers, are 
essential for stable operation of BASIS, and are used to create, maintain, and modify its 
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analog control levels (the control levels of the comparators in the Threshold Domain Filter). 
For example, a baseline created by an AARF operating as a median filter (i.e., q = 1/2) 
will not significantly change its value unless the photoelectron rate exceeds about half of the 
saturation rate i^max- On the other hand, this baseline will track the changes in the noise 

5 level, providing an effective separation between noise and the photosignal. 

When the photoelectron rate exceeds the saturation rate i2max> the output of the AARF 
itself will well represent the central tendency of the photosignal, and thus will be proportional 
to the incident photon rate. In the 'transitional' region (around i2max), the output of the 
BASIS can be constructed as a weighted sum of the outputs of AARF and ACM. Thus the 

10 total output of BASIS can be constructed as a combination of the outputs of AARF, ACM, 
and SRM, and calibrated to be proportional to the incident photon rate. 
Analog Counting Module (ACM) This module produces a continuous output, R(t), equal to 
the rate of upward zero crossings of the difference, x(t) — rD q (t; T), in the time window, w(t), 
given by 



15 R(t) = w(t) * K(t) = w* 



(43) 



where 7£(i) denotes the instantaneous crossing rate (Nikitin et al., 2003). The value of the 
parameter r generally depends on the distribution of the photosensor's noise in relation to 
the single photoelectron distribution of the photosensor, and can normally be found either 
theoretically or empirically based on the required specifications. This parameter affects the 

20 ratio of the false positive (noise) and the false negative counts (missed photoelectrons) and 
allows us to achieve a desired compromise between robustness and selectivity. In the subse- 
quent simulated example (see figure 25), the quantile parameter q = 1/2 (AARF operating 
as a median filter), and r = 6. An attractive choice for the time averaging filter w(t) is a 
sequence of 3 RC-integrators with identical time constants r — T/6, which will provide us with 

25 rate measurements corresponding to the time averaging with a rectangular moving window of 
width T (Nikitin et ai., 2003). 

The main advantage of the analog counting represented by equation (43) is a complete 
absence of dead time effects (see Nikitin et al. 2003). In addition, the baseline created by 
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an AARF will not be significantly affected by the photoelectron rates below approximately 
(1 — q) Rmax (half of the photosensor saturation rate for a median filter). Thus the maximum 
measured rate is limited only by the single electron response of a photosensor. This is at least 
two orders of magnitude higher than the current state of the art photon counting systems. For 
5 example, in the simulation presented in figure 25, the FWHM of the single electron response 
is about 1 ns, and the saturation rate of photon counting is about 3 x 10 8 s" 1 . Since the 
signal-to-noise ratio is proportional to the square root of the rate, the 20 dB increase in the 
photon counting rate translates into a tenfold increase either in the sensitivity or the speed of 
detection. 

10 Saturation Rate Monitor (SRM) The SRM produces a continuous output -Rmax(0 equal to 
the rate of upward zero crossings of the difference x(t) — D 1 /2(t] T) in the time window w(t) 1 



5 «(-«,) 



(44) 



As was theoretically derived by Nikitin et al. (1998), iimax is approximately equal to the 

maximum rate of upward (or downward) crossings of any constant threshold by the pho- 
15 tosensor signal x(t). When the photoelectron rate Xphe of a photosensor is much smaller 

than .ftmax) the pileup effects are small, and the photosensor is in a photon counting mode. 

When Xphe flmax, the photosensor is in a current mode. 

Thus monitoring i?max allows us to automatically handle the transition between the two 

modes. The horizontal gray line in panel I of figure 25 shows the measured Rmax as a function 
20 of the photoelectron rate Xphe- The measured saturation rate is also shown by the horizontal 

thin solid lines in the lower half of panel III of figure 25. 

Integrated Output Module As shown in figure 24, the output module of BASIS combines 
the outputs of the ACM, SRM, and AARF into a single continuous output Ro U t(t). The 
magnitude of Rout(t)> for both the photon counting and the current modes, is proportional to 
25 the rate of incident photons. This signal can be, for example, used directly in analog or digital 
measuring and/or control systems, differentiated (thus producing continuous time derivative 
of the incident photon rate), or digitally sampled for the subsequent transmission and/or 
storage. For the simulated example shown in figure 25, Ro Ut (t) was chosen as the following 
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combination of D q (t;T), R(t), and #maxC0 : 

flout(t) = R(t) + PD q (t; T) T* D \pD q (t\ T) - 7 i2max(i)] , (45) 

where (3 is a calibration constant, AD = aRm^^ a being a small number (of order 10" 1 ), and 
7 ~ 1/2 is a quantile constant. The Integrated Output Module thus includes the 'transitional' 
5 region between the photon counting and the current modes (shaded in gray in figure 25), 
currently unavailable, into a normal operational range of a photosensor, extending by ~20 dB 
the photosensor's contiguous dynamic range. 

Simulated examples of light measurements conducted by PMT with BASIS unit Figure 25 
provides a simulated example of the performance of BASIS used with a PMT. In the sim- 

10 ulation, a fast PMT was used (the FWHM of the single electron response is about 1 ns), and 
the noise rate was chosen to be high (order of magnitude higher than the PMT saturation 
rate). Panel I of the figure shows the output of BASIS (i?o U t, thick solid black line) as a 
function of the photoelectron rate Ap h E, along with the outputs of the Saturation Rate Moni- 
tor (iCax, solid gray line), Rank Filtering Module (£>i/2> dashed line), and Analog Counting 

15 Module (i?, thin solid black line). Panel II shows (by gray lines) 1 /xs snapshots of the PMT 
signal for the photoelectron rates much smaller (left), approximately equal (middle), and much 
higher (right) than the saturation rate of the PMT. This panel also shows (by the black lines) 
the respective outputs of the Rank Filtering Module Di/ 2 (t), and the baseline levels rDi/ 2 (t) 
used in the Analog Counting Module (r = 6 in the simulation). The instantaneous crossing 

20 rates H(t) are also shown (top), and the time constant T of AARF is indicated in the lower left 
corner of the panel. Panel III illustrates the relation between the noise and photosignal used 
in the simulation by depicting the accumulated amplitude and counting distributions of the 
PMT signal. These distributions are shown for three chosen photoelectron rates A PhE , in their 
relation to the outputs of the Rank Filtering Module (Di/ 2 ) and the Saturation Rate Monitor 

25 (Rmax)- The resolution AD of the acquisition system used for measuring the distributions is 
indicated in the panel. 

Figure 26 provides a simulated example of a modification of BASIS designed for detection 
of fast changes in a light level. The light signal corresponding to this model can be, for 
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example, an intensity modulated light signal passing through a fiber, or fluorescence of dye 
excited by an action potential wave propagating through a biological tissue. The gray line in 
the lower panel of the figure shows the time- varying light signal (square pulses). The higher 
light level corresponds to the photoelectron rate of about 2 x 10 9 photoelectrons per second. 
The width (FWHM) of the single electron response of the photosensor is about 1 ns, and the 
resulting photosensor electrical signal x(t) is shown by the gray line in the middle panel. 

As can be seen in the figure, the low signal-to-noise ratio makes fast and accurate deduction 
of the underlying light signal difficult. The circuit shown at the top of figure 26, however, 
allows reliable timing of the onsets and offsets of the light pulses with better than 10 ns 
accuracy. The output of the circuit D q{ (t\T{) is shown by the black line in the lower panel 
of the figure. In the example, the quantile parameters of the rank filters are g& = 1/4 and 
q f = 3/4, and the baseline factor is r = 1.5. The parameter r allows us to adjust the circuit 
for optimal performance based on the difference between the low and high light levels. 
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5 Generation of monoenergetic Poissonian pulse trains 

As another illustration of the current invention, consider a technique and a circuit for genera- 
tion of monoenergetic Poissonian pulse trains with adjustable rate and amplitude. Generators 
of such pulse trains can be used, for example, in testbench development and hardware proto- 

5 typing of instrumentation for nuclear radiation measurements. 

Idealized model of a Poisson pulse train generator An idealized process producing a monoen- 
ergetic Poissonian pulse train can be implemented as schematically shown in figure 27. Con- 
sider a stationary random pulse train J2% %i S(t — where X{ are the amplitudes of the pulses 
with the arrival times U. This pulse train is filtered by a linear time filter with a continuous 

10 impulse response w& T (t), where At is the characteristic response time of the filter. An exam- 
ple of such a response would be the bipolar pulse WA T (t) = [i/Ar 2 ~ t 2 /(2Ar 3 )] e~ </Ar 0(i), 
where 6 is the Heaviside unit step function. The output x(t) of the linear time filter can be 
written as 

x(t) =Y^Xi w Ar (t - U) , (46) 

i 

15 and is a continuous signal. The instantaneous rate of upward crossings Nikitin et al. (2003) 
of a threshold D by this signal can be written as 

= E*(*-*j). (47) 

+ j 

where tj are the instances of the crossings (that is, x(tj) = D and x(tj) > 0). As was 
discussed in Nikitin (1998), the pulse train given by equation (47) is an approximately Pois- 
20 sonian train affected by a non-extended dead time of order Thus, when the average 

rate R(D) — {1Z(D,t)) T is much smaller than the saturation rate iimax, TZ(D,t) will provide 
a good approximation for a monoenergetic Poissonian pulse train of the average rate R(D). 

When either Mi = 0 or Wig — 0, then, as was shown in Nikitin et al. (1998), the average 
rate of the upward crossings of a threshold D by the signal x(t) can be expressed as 



K{D,t) = 



-8[x{t)-D] 
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R(D) = i^ax exp j-± (£)*] , (48) 

and thus the rate of the generated pulse train can be adjusted by an appropriate choice of the 
threshold value D. 

Practical implementation of a Poisson pulse train generator The idealized process described 
above is not well suited for a practical generation of a Poissonian pulse train, since, as can be 
seen from equation (48), at high values of the threshold D the rate of the generated train is 
highly sensitive to the changes in D. To reduce this sensitivity, one can pass the signal x(t) 
through a nonlinear amplifier, e.g., an antilogarithmic amplifier as shown in figure 27, thus 
transforming x(t) into the signal y(t) = exp [x(t)/a]. Then the average rate of the upward 
crossings of a threshold D by the signal y(t) can be written as 

^^i^axexpj-ijln^)] 2 ! , (49) 

which is much less sensitive to the relative errors in D. 

Figure 28 illustrates a simulated performance of an idealized monoenergetic Poisson pulse 
generator shown in figure 27. The upper panel of figure 28 shows the output pulse rates as a 
function of threshold, and the lower panels show the distributions of the pulses' interarrival 
times for the generator set at two different threshold values. In the figure, the black solid 
lines show the theoretical curves, and the gray solid lines show the respective results of the 
simulations. 
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Articles of Manufacture 

Various embodiments of the invention may include hardware, firmware, and software em- 
bodiments, that is, may be wholly constructed with hardware components, programmed into 
firmware, or be implemented in the form of a computer program code. 
5 Still further, the invention disclosed herein may take the form of an article of manufacture. 

For example, such an article of manufacture can be a computer-usable medium containing a 
computer-readable code which causes a computer to execute the inventive method. 
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Regarding the invention being thus described, it will be obvious that the same may be 
varied in many ways. Such variations are not to be regarded as a departure from the spirit 
and scope of the invention, and all such modifications as would be obvious to one skilled in 
the art are intended to be included within the scope of the following claims. 

We claim: 

1. Method for processing an input signal consisting of a Threshold Domain Filtering 
step utilizing a threshold domain and at least one of the following additional steps: 

(a) Multimodal Pulse Shaping; 

(b) Analog Rank Filtering; and 

(c) Analog Counting. 

2. Method for signal processing of claim 1 wherein said Threshold Domain Filtering 
step is used for extraction of one or more features of interest from said input signal 
to produce one or more output signals. 

3. Method for signal processing of claim 1 wherein said threshold domain is defined by 
a Control Level Signal consisting of one or more components. 

4. Method for signal processing of claim 1 wherein said Threshold Domain Filtering step 
outputs a first output value when said input signal is within said threshold domain 
and a second output value when said input signal is outside of said threshold domain. 

5. Method for signal processing of claim 1 wherein said Threshold Domain Filtering step 
outputs a first output value when said input signal is within said threshold domain, a 
second output value when said input signal is outside of said threshold domain, and 
one or more values when said input signal is inside a finite interval at a boundary of 
said threshold domain. 
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6. Method for signal processing of claim 1 wherein at least one of said additional steps 
is Multimodal Pulse Shaping and wherein said input signal has one or more com- 
ponents and wherein said Multimodal Pulse Shaping transforms at least one of said 
components of said input signal into a plurality of signal components where at least 
one said component represents a derivative of another said component. 

7. Method for signal processing of claim 6 wherein said derivative is a time-derivative. 

8. Method for signal processing of claim 6 wherein said plurality of signal components 
is representative of one or more features of interest of said input signal. 

9. Method for signal processing of claim 3 wherein one or more components of said 
Control Level Signal are outputs of a step of Analog Rank Filtering. 

10. Method for signal processing of claim 3 wherein one or more components of said 
Control Level Signal are results of a step of Analog Rank Filtering of one or more 
output components of a step of Multimodal Pulse Shaping of said input signal. 

11. Method for signal processing of claim 1 wherein at least one of said additional steps is 
Analog Counting, wherein said Threshold Domain Filtering step has an output with 
at least one component and where said Analog Counting comprises the steps of: 

(a) time-differentiating said component to produce a time-derivative, and 

(b) rectifying said time-derivative to produce an instantaneous count rate. 

12. Method for signal processing of claim 11 wherein said Analog Counting further com- 
prises the step of: 

(a) time-integrating said instantaneous count rate to produce a count rate in a mov- 
ing time window. 
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13. Method for Adaptive Analog Rank Filtering operable to transform an input scalar 
variable into an output Rank Filtered Variable comprising steps of: 

(a) forming a plurality of outputs of delayed comparators by passing said input 
variable and a plurality of feedbacks of Offset Rank Filtered Variables through 
a respective plurality of said delayed comparators; 

(b) forming a first difference which is a weighted difference of said outputs of the 
delayed comparators; 

(c) forming a Density Function by time averaging said first difference; 

(d) forming a second difference which is a weighted difference of said feedbacks of 
the Offset Rank Filtered Variables; 

(e) forming a plurality of differences between the outputs of the delayed compara- 
tors and the respective Offset Quantile Parameters of said Offset Rank Filtered 
Variables; 

(f) forming a plurality of time derivatives of said Offset Rank Filtered Variables by 
multiplying each of said plurality of differences by a ratio of said second difference 
and said Density Function; 

(g) producing the plurality of the Offset Rank Filtered Variables by time-integrating 
said plurality of time derivatives; and 

(h) producing said Rank Filtered Variable as a weighted average of said Offset Rank 
Filtered Variables. 
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14. Method for Adaptive Analog Rank Selecting operable to transform a plurality of input 
variables into an output Rank Selected Variable comprising steps of: 

(a) forming a plurality of outputs of averaging comparators by passing said plurality 
of input variables and each component of a plurality of feedbacks of Offset Rank 
Selected Variables through a respective plurality of said averaging comparators; 

(b) forming a first difference which is a weighted difference of said outputs of the 
averaging comparators; 

(c) forming a Density Function by time averaging said first difference; 

(d) forming a second difference which is a weighted difference of said feedbacks of 
the Offset Rank Selected Variables; 

(e) forming a plurality of differences between the outputs of the averaging compara- 
tors and the respective Offset Quantile Parameters of said Offset Rank Selected 
Variables; 

(f) forming a plurality of time derivatives of said Offset Rank Selected Variables by 
multiplying each of said plurality of differences by a ratio of said second difference 
and said Density Function; 

(g) producing the plurality of the Offset Rank Selected Variables by time-integrating 
said plurality of the time derivatives; and 

(h) producing said Rank Selected Variable as a weighted average of said Offset Rank 
Selected Variables. 
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15. Method for signal processing of claim 3 wherein one or more components of said 
Control Level Signal are linear combinations of outputs of a step of Adaptive Analog 
Rank Filtering of one or more components of said input signal. 

16. Method for signal processing of claim 3 wherein one or more components of said 
Control Level Signal are linear combinations of outputs of a step of Adaptive Analog 
Rank Filtering of one or more output components of a step of Multimodal Pulse 
Shaping of said input signal. 

17. Method for signal processing of claim 3 wherein one or more components of said 
Control Level Signal are linear combinations of outputs of a step of Adaptive Analog 
Rank Selecting of a plurality of components of said input signal. 

18. Method for signal processing of claim 3 wherein one or more components of said 
Control Level Signal are linear combinations of outputs of a step of Adaptive Analog 
Rank Selecting of a plurality of output components of a step of Multimodal Pulse 
Shaping of said input signal. 
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